Finite-size correction in many-body electronic structure calculations 
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Finite-size (FS) effects are a major source of error in many-body (MB) electronic structure cal- 
culations of extended systems. A method is presented to correct for such errors. We show that 
MB FS effects can be effectively included in a modified local density approximation calculation. A 
parametrization for the FS exchange-correlation functional is obtained. The method is simple and 
gives post-processing corrections that can be applied to any MB results. Applications to a model 
insulator (P2 in a supercell), to semiconducting Si, and to metallic Na show that the method delivers 
greatly improved FS corrections. 
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Realistic many-body (MB) calculations for extended 
systems are needed to accurately treat systems where 
the otherwise successful density functional theory (DFT) 
approach fails. Examples range from strongly correlated 
materials, such as high-temperature superconductors, to 
systems with moderate correlation, for instance where 
accurate treatments of bond-stretching or bond-breaking 
are required. DFT or Hartree Fock (HF), which are ef- 
fectively independent-particle methods, routinely exploit 
Bloch's theorem in calculations for extended systems. In 
crystalline materials, the cost of the calculations depends 
only on the number of atoms in the periodic cell, and the 
macroscopic limit is achieved by a quadrature in the Bril- 
louin zone, using a finite number of k-points. MB meth- 
ods, by contrast, cannot avail themselves of this simpli- 
fication. Instead calculations must be performed using 
increasingly larger simulation cells (supercells). Because 
the Coulomb interactions are long-ranged, finite-size (FS) 
effects tend to persist to large system sizes, making reli- 
able extrapolations impractical. The resulting FS errors 
in state-of-the-art MB quantum simulations often can be 
more significant than statistical and other systematic er- 
rors. Reducing FS errors is thus a key to broader appli- 
cations of MB methods in real materials, and the subject 
has drawn considerable attention [l], Q • 

In this paper, we introduce an external correction 
method, which is designed to approximately include FS 
corrections in modified DFT calculations with finite-size 
functionals. The method is simple, and provides post- 
processing corrections applicable to any previously ob- 
tained MB results. Conceptually, it gives a consistent 
framework for relating FS effects in MB and DFT calcu- 
lations, which is important if the two methods are to be 
seamlessly interfaced to bridge length scales. The cor- 
rection method is applied to a model insulator (P2 in a 
supercell), to semiconducting bulk Si, and to Na metal. 
We find that it consistently removes most of the FS er- 
rors, leading to rapid convergence of the MB results to 
the infinite system. 



cell as (Rydberg atomic units are used throughout): 

N M 

H = E v ? + E *w + E ^ FS (i r * - ^i) . 

i— 1 i—1 i<j 

where the ionic potential on i can be local or non-local, 
and Yi is an electron position. The Coulomb interaction 
V between electrons depends on the supercell size and 
shape, due to modification by the periodic boundary con- 
ditions (PBC) Q. A FS correction is often applied to the 
MB results from parallel DFT or HF calculations. The 
corresponding DFT, as usually formulated, introduces a 
fictitious mean-field A/"-electron system 0, [a] : 



DFT 



-V 2 + V ion + V H {r) + V™{r) , 



(2) 



where the Hartree and exchange-correlation (XC) poten- 
tials depend self-consistently on the electronic density 
71 (r). In the non-spin-polarized local density approxima- 
tion (LDA), for example: V£{r) = 5{n{r) ef c {n))/8n{v), 
where e%? c (n) is typically obtained from quantum Monte 
Carlo (QMC) results on the homogeneous electron gas 
(jellium), extrapolated to infinite size [1, 0, @]- 

Residual errors after DFT FS correction are still found 
to be large, however, and the equations above illustrate 
why. The jellium QMC results, which determine e^(n), 
have been extrapolated to infinite supercell size for each 
density. This is the correct choice for standard LDA ap- 
plications, where Bloch's theorem will be used to reach 
the infinite limit. It is not ideal, however, if the LDA 
is expected to provide FS corrections. Only one-body 
FS corrections (kinetic, Hartree, etc), which arise from 
incomplete k-point integration, are included, while two- 
body FS corrections [l| are missing. If parallel HF cal- 
culations are used instead, exact FS exchange V x — > 
is included, but V c is zero. 

Our approach is to construct an LDA with FS XC in 
Eq. @ . If the supercell of Eq. |T]) is cubic (for simplicity) , 
the XC energy is e£f (n) = e x (r s ,L) + e c (r s ,L), where r s 
specifies the density via 47rrf/3 = 1/n and L denotes 
the linear size of the supercell. To obtain e^f(r s ,L), we 
We write the A/"-elcctron MB Hamiltonian in a super- use unpolarized jellium systems in the same supercell, in 
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FIG. 1: (Color online) Calculated and parametrized jellium 
exchange and correlation energies per electron vs. r S) for a 
range of supercell sizes L (in Bohr). Top panel: exchange 
energies, with solid lines given by the fit in Eq. ((3}, and open 
symbols by k-point averaged calculations. Bottom panel: cor- 
relation energies, with solid lines given by the fit in Eq. ((6}, 
open symbols by AF QMC calculations, and small filled sym- 
bols by the large- TV asymptotic expression in Eq. <(Sj . 



which the number of electrons N (distinct from N) is a 
variable given by the ratio r s /L: N = (3/47r)(L/r s ) 3 . 
We parametrize the HF exchange energy in jellium by: 



e x (r s ,L) = 



Hi, 



if r s < 7; 
other. 



(3) 



The term with clq ~ — 0.916 Ry gives the usual infinite- 
size limit. A 1/L canceling term, which arises from the 
self- interaction of an electron with its periodic images [|[ , 
has been implicitly included. The leading FS dependence 
is then 1/L 2 0] . The form of the remaining terms is moti- 
vated by the exact scaling relation: e x (r s , L) = e x (N)/L. 
To obtain e x , we calculate e x (N) for a range of N, each 
by averaging over about 20 k-points. The results are fit- 
ted to give ai and 0.2- As illustrated in Fig.[Tl the quality 
of the fit is excellent. The behavior of e x at large r s re- 
quires special handling for finite L. At 7 = r s (N = 2), 
there is only one electron of each spin in the supercell, so 
e x is just the self-interaction term. Beyond 7, e x is forced 
to go to zero as 1/Vf, reflecting the self- interaction of a 
'fractional' electron. The coefficient 03 is chosen to make 
the exchange potential V x s continuous at 7 ■ From 
Fig. [1] the magnitude of the discontinuity at r s = 7 is 
seen to decrease with increasing L, as expected. All pa- 
rameters are listed in Table |TJ 

The correlation energy in jellium is the difference be- 
tween the MB and HF energies (per electron): 



TABLE I: Parameters (in Ry atomic units) in the FS XC 
functionals. 
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[Eq. ©] 


-2.2037 0.4710 


-0.0150 
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[Eq. ©] 


0.1182 1.1656 


-5.2884 


-1.1233 



where the jellium non-interacting kinetic energy obeys 
the scaling relation t(r s ,L) = t(N)/L 2 . We calculate 
t(N) in the same way as e x (r s ,L), but averaging over 
more k-points to ensure convergence. 

We next derive the MB energy £ . Ceperley and Alder 
obtained jellium QMC energies for various values of 
N and provided the following fit: 

S(r s ,L)=e co (r s ) + b 1 (r s )At N + b 2 (r s )/N, (5) 

where N uniquely determines L, and Atjy = t(r s ,L) — 
t(r s , 00) is the FS error in the free-electron kinetic energy. 
The infinite-size limit, £°°, was extrapolated from Eq. ([5]) 
and it is the basis for V™ in Eq. The b parameters 
were given for several r s values, which we fit to get the 
functions b\ (r s ) and 6 2 (r s ). With these and At N , we can 
now calculate £(r Sl L) for any r s and L, which is accurate 
for large N. 

For small N, namely large r s in a finite supercell, 
Eq. ([5]) does not apply. This is easy to see from the 
1/N term which, at sufficiently large r s , causes e c to di- 
verge. To guide the analysis in this region, we use the 
plane- wave auxiliary-field (AF) QMC method [H Q3] to 
directly calculate £ for FS jellium systems. At small r s , 
the correlation energy obtained is in excellent agreement 
with that derived from Eq. ((5]), as shown in Fig. Q] At 
large r s , e c from Eq. ([5|) falls below the AF QMC value 
[b2{r s ) is negative], as the latter goes to zero monoton- 
ically. The value of r s where the two begin deviating 
depends on L, since it is determined by N. 

We thus parametrize the correlation energy by 



0, 



r s < 7h\ 
lh <r s < 7; ; 
other. 



e c (r Sl L) = £{r s .L) - t(r s ,L) - e x (r s ,L), 



(4) 



(6) 

The correlation functional has been divided into high, 
intermediate, and low density regions. The boundaries 
are defined by 7/j = r s (N = 12) and 7/ = r s (N = 
1/2), which are guided by the discussion in the pre- 
vious paragraph and the quality of the fits described 
below, but are otherwise arbitrary. At high densities, 
the infinite-size limit is given by e£°(r s ) (the Perdew- 
Zunger parametrization [8| is used here), and the lead- 
ing FS term exactly cancels that in e x (r s ,L), to ensure 
that e xc (r s , L) correctly scales as 0(1/ L 3 ). The function 

g(r s ) = g\r s ln(r s ) +gzr s + g^V 2 + gtr 2 is obtained from 
a fit to e c (r s ,L) from Eqs. (g]) and (0. (The fits are il- 
lustrated in Fig. [T] and parameters are given in Table [T]) 
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At intermediate densities, the function f(r s ) is given by 
a cubic polynomial and is completely determined by the 
requirement that e c and its derivative be continuous at 
r s = Ih an d t s =7;. As Fig. [1] shows, the parametriza- 
tion in Eq. ^ closely reproduces our AF QMC data at 
low densities for all cell sizes. 

Post-processing FS corrections are now easily gener- 
ated for any MB calculation. The DFT FS results, using 
e£f from Eqs. ((3]) and (|6|), can be obtained from stan- 
dard DFT computer codes with only minor modifica- 
tions. If E FS (L) is the energy from DFT FS and E(L) 
from standard DFT (i.e., DFT°°), the energy correction 
is ADFT FS = E(oo) - E FS (L), where E(oo) is obtained 
by k-point integration. The correction can alternatively 
be expressed as the sum of ADFT 1B = E{po)-E{L) and 
ADFT 2B = E(L) - E FS (L). The one-body (IB) correc- 
tion is the usual ADFT°°, while the two-body (2B) part 
captures the FS effects that arise from the modification 
of V xc due to supercell PBC. 

The present correction scheme is exact for homoge- 
neous systems. Our first application of the method is 
to a model system in the opposite limit. We consider a 
"molecular solid" with P2 in a periodic supercell, treated 
by the plane- wave AF QMC method Because of 

the low-density "vapor" region and the variation in den- 
sity, the system provides a challenging test for the cor- 
rection method. A norm-conserving Klcinman-Bylandcr 
13 1 separable non-local LDA pseudopotential is used 
14j . Total energy calculations were performed at the 
equilibrium bondlength of 3.578 Bohr, for cubic super- 
cells of size L = 7—18 Bohr, all with the T-point (k = 0). 
Figure [2] shows the results from AF QMC and LDA us- 
ing both DFT°° and DFT FS [15]. The uncorrected QMC 
result has large FS errors and, at L = 18, is still ~ 0.3 eV 
away from the infinite-size value. Corrected with DFT°°, 
the FS error is somewhat reduced at intermediate L, but 
is unchanged for larger L where the 2B effects dominate. 
With the new method, the corrected energy shows excel- 
lent convergence across the range, reaching the asymp- 
totic value (within statistical errors) by L ~ 12. 

The second application is for fee bulk silicon, using 
non-cubic supercells (nxnxn the size of the primitive 
fee cell). The raw MB energies in Fig. [3] are taken from 
diffusion Monte Carlo (DMC) calculations Q. In the FS 
corrections for the fee supercells, we use e F f from Eqs. ^ 
and ([6]), with an effective L equal to the size of a cubic 
supercell of the same volume. The pseudopotential used 
is also different from that in the DMC calculations. We 
checked multiple pseudopotentials to ensure that the FS 
corrections are independent of the choice of pseudopoten- 
tial. The DMC calculations were done with the k = L 
point. The usual DFT correction is in the wrong direc- 
tion in this case, thereby increasing the FS error. The 
new method removes most of the error, despite the non- 
optimal e F f . The inset in Fig. [3] shows ADFT 2B cal- 
culated as above for fee, compared with that for cubic 
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FIG. 2: (Color online) P2 PBC total energy convergence 
vs. supercell size. Standard DFT°° FS effect is different (too 
small) from that of MB AF QMC. DFT FS parallels the MB 
calculation and leads to much more rapid convergence. The 
inset focuses on larger L and shows the raw and corrected AF 
QMC results plotted as a function of 1/L 3 . 
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FIG. 3: (Color online) Total energy per atom (in eV) of bulk 
Si in fee supercells of size nxnxn (2n 3 atoms). DMC energies 
in the main graph are from Kent et. al. [l[ (shifted here rela- 
tive to the extrapolated infinite-size limit). FS corrections are 
shown from both standard DFT and the present method. The 
inset shows the calculated two-body correction as a function 
of inverse the number of atoms, for fee and cubic supercells. 



supercells. Both are seen to fall on an essentially smooth 
and linear curve. This weak shape dependence of e F f is 
encouraging, suggesting that additional FS MB jcllium 
calculations can be avoided in some non-cubic supercells. 

The final application is for metallic bec bulk Na. While 
in insulators a single k-point is often adequate, metals 
present additional difficulties. We do multiple MB cal- 
culations with random kjxrints (e.g. 50 for 16 atoms) 
and average the results [16[. The plane- wave AF QMC 
method (ill . Il2| was used, in which any k-point can be 
included by a simple modification to the one-particle ba- 



TABLE II: Calculated cohesive energy (in eV) of bcc solid 
sodium vs. experiment. AF QMC results from supercells with 
2, 16 and 54 atoms are shown, together with FS-corrected 
results. A zero-point energy of 0.0145 eV/atom is included. 
Experimental value was taken from Ref. [171 ] . 







corrected 




raw 


w/ 1-body w/ full FS 


2 


2.050(35) 


2.141(2) 1.124(2) 


16 


1.264(14) 


1.287(4) 1.135(4) 


54 


1.184(9) 


1.189(10) 1.143(10) 


expt 




1.129(6) 



sis. Although our pseudopotential has a Ne-core, DFT 
tests with various pseudopotentials verified that it is suf- 
ficient for the cohesive energy (consistent with Ref. [l7l]). 
but the frozen semi-core introduces systematic biases in 
the lattice constant and bulk modulus. The calculated 
cohesive energies are given in Table HU The FS-corrected 
cohesive energies for 16 and 54 atoms are consistent, and 
in better agreement with experiment than the previous 
best DMC results [13] of 1.0221(3) eV (with a core polar- 
ization potential) and 0.9910(5) eV (without). The cal- 
culated equation of state is shown in Fig. [4] We see 
that, with the new FS corrections and k-point sampling, 
the calculations have better convergence than previously 
reachable with an order of magnitude larger system sizes 
17| ■ Both the lattice constant and the bulk modulus were 
modified by the FS corrections. As the bottom panel 
demonstrates, FS effects always cause a systematic error 
in the lattice constant in uncorrected MB calculations. 

These tests show that our DFT FS correction method 
works well in a variety of systems. This is perhaps not 
surprising, given the often near-sighted nature of the XC 
function. For the method to be effective, DFT needs to 
provide a good approximation in capturing the differ- 
ence between the systems with interaction V FS and V, 
which is not the same as requiring DFT to work well in 
either system (assuming L greater than the size of the 
XC hole). We have presented an XC functional which 
delivers high accuracy across several different materials. 
Previous attempts at FS correction have focused on es- 
timating the errors internally within the MB simulation 
P, [3] ■ Our approach is an external method which is sim- 
ple and can provide post-processing FS correction to any 
MB electronic structure calculations. The method can 
be generalized, e.g., to spin-polarized systems and other 
supercell shapes, and the FS functional could be further 
improved, e.g., by exact exchange. 

We thank E. J. Walter for help with pseudopoten- 
tials, W. Purwanto for help with computing issues, and 
P. Kent for sending us the numerical data from Ref. []]]. 
This work is supported by ONR (N000140510055), NSF 
(DMR-0535529), and ARO (48752PH) grants. Comput- 
ing was done on NERSC and CPD computers. 
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FIG. 4: (Color online) Top: Equation of state for bcc Na. The 
AF QMC energy is shown vs. lattice constant a for 16- and 54- 
atom supercells, with one-body (IB) and full FS corrections. 
With full correction (solid lines), 16- and 54-atom results are 
almost indistinguishable. QMC statistical errors include that 
of k-point averaging. The vertical arrows indicate the calcu- 
lated equilibrium a. Bottom: persistence of the two-body FS 
error, with finite slopes (number of atoms indicated). 
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